Gravitational wave snapshots of generic extreme mass ratio inspirals 



Steve Drasco 

Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109 and 
Center for Radiophysics and Space Research, Cornell University, Ithaca, NY 14853 

Scott A. Hughes 
Department of Physics and MIT Kavli Institute, 
MIT, 77 Massachusetts Ave., Cambridge, MA 02139 

Using black hole perturbation theory, we calculate the gravitational waves produced by test 
particles moving on bound geodesic orbits about rotating black holes. The orbits we consider 
are generic — simultaneously eccentric and inclined. The waves can be described as having radial, 
polar, and azimuthal "voices", each of which can be made to dominate by varying eccentricity 
and inclination. Although each voice is generally apparent in the waveform, the radial voice is 
prone to overpowering the others. We also compute the radiative fluxes of energy and axial angular 
momentum at infinity and through the event horizon. These fluxes, coupled to a prescription for 
the radiative evolution of the Carter constant, will be used in future work to adiabatically evolve 
through a sequence of generic orbits. This will enable the calculation of inspiral waveforms that, 
while lacking certain important features, will approximate those expected from astrophysical extreme 
mass ratio captures sufficiently well to aid development of measurement algorithms on a relatively 
short timescale. 
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I. INTRODUCTION 

One of the leading problems in modern general rela- 
tivity research is the analysis of binary systems. This 
problem is quite nasty in general, especially when the 
members of the binary are strongly interacting and 
gravitational-wave (GW) emission is significant. It sim- 
plifies considerably in the extreme mass ratio limit, in 
which one member of the binary is much smaller than 
the other (taken to be a black hole). This is largely be- 
cause the small mass ratio makes it possible to bring 
perturbative techniques to bear. 

This paper reports the development of a key piece of in- 
frastructure needed for an approximate approach to this 
problem: the ability to accurately calculate the instanta- 
neous GW emission from arbitary bound geodesic orbits 
of rotating Kerr black holes. Our work generalizes previ- 
ous calculations which specialized the orbit to be either 
circular (of constant Boyer-Lindquist coordinate radius 
0) or equatorial (confined to the black hole's plane of 
reflection symmetry |2j). 

By knowing the GWs, we also calculate the flux of 
energy E and axial angular momentum L z that the radi- 
ation carries to infinity and down the black hole's event 
horizon. Global conservation allows us to evolve the en- 
ergy and angular momentum associated with the orbit, 
fixing the evolution of two of the three "constants" which 
characterize black hole orbits (up to initial conditions). 
The remaining constant, Q, is known as the Carter con- 
stant 1 0, 3 . In the limit of zero black hole spin, the hole 



1 The constant K = Q + (L z — aE) 2 is often used as well. 



is spherically symmetric, geodesic orbits conserve total 
angular momentum L, and Carter's constant is given by 
Q = L 2 — L 2 Z ; though not strictly accurate, the intu- 
itive interpretation "Q = the rest of the angular mo- 
mentum" is useful even for non-zero black hole spin. It 
has recently become clear that, when the system evolves 
"slowly" (quantified below), it is possible to fix the evo- 
lution of Q using techniques not too different from those 
used to evolve E and L z USUI!. 

For such slowly evolving systems, we can then evolve 
all of the orbital constants; by doing so, we can build rea- 
sonably realistic inspirals corresponding to generic initial 
conditions, as well as the associated waveforms. The 
qualifier "reasonably realistic" refers to the fact that 
this procedure neglects, by construction, the influence of 
"conservative" pieces of the small body's self interaction; 
only the "dissipative" influence can be analyzed in this 
way. As we discuss further below, neglect of the conser- 
vative influence will limit the reliability of our inspirals; 
for the applications we have in mind, this limitation is 
a reasonable trade off in order to build waveforms on a 
rapid timescale. 

We now briefly summarize our major motivations for 
analyzing this problem, sketch the techniques that we 
have developed and our major results, and discuss the 
utility and limitations of our results. 



A. Astrophysical context 

Though an idealization of a much more difficult general 
problem, the extreme mass ratio limit in fact corresponds 
precisely to an important subset of astrophysical bina- 
ries: stellar mass compact objects captured onto highly 
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relativistic orbits of massive black holes in galaxy cores. 
Such binaries are expected to be created by multibody 
scattering processes in the stellar cluster that surrounds 
these central black holes. 

It has now become clear that the central bulges of es- 
sentially all galaxies contain a massive black hole with 
properties strongly correlated to that of the bulge P.ITo|. 
It also appears likely that mass segregation (the sinking of 
massive stellar objects to the bottom of a gravitational 
potential due to equipartition of kinetic energy is 
likely to drive the largest stellar mass compact objects 
into the vicinity of the massive black hole, increasing 
the probability of forming a capture binary — see, for 
example, Ref. [IsJ f° r discussion pointing to the possi- 
bility of a "minicluster" of stellar mass black holes in 
the center of our galaxy, and Ref. [l3| for evidence of an 
overabundance of black holes in binaries in the galaxy's 
inner parsec. Mass segregation is also seen numerically in 
galactic center N-body simulations by Freitag 14]. The 
abundance of galactic black holes plus mass segregation 
suggest that the capture and GW driven inspiral of stel- 
lar mass compact objects into massive black holes may 
be a relatively common phenomenon in the universe. 

The GWs produced by these extreme mass ratio inspi- 
rals (EMRIs) are thus ideal targets for low-frequency GW 
detectors, particularly the space-based LISA antenna, 
currently under development as a joint NASA-ESA mis- 
sion (la . Ha ]. By combining estimates for the rate at 
which such captures are likely to occur (per galaxy) with 
the projected sensitivity of LISA (which determines the 
distance to which sources can be measured, and hence 
sets the number of galaxies that are relevant), one finds 
that dozens to thousands of EMRI events may be mea- 
sured over LISA's mission lifetime [Tt| . 

Because the small body only slightly perturbs the 
spacetime of the large black hole, the GWs are largely 
set by — and therefore encode — the nature of the black 
hole's spacetime. A typical measured EMRI event is 
likely to last for ~ 10 4 — 10 5 orbits; by coherently tracking 
the waveform's phase over these orbits, LISA should de- 
termine parameters characterizing the binary with great 
precision. Black hole masses and spins should be mea- 
sured with fraction of a percent accuracy or better |l8| ; it 
should even be possible to check whether the spacetime 
has the "shape" (higher multipolar structu re) requ ired 
by the no-hair theorems of general relativity (lil 120. 12H] . 
EMRI measurements are expected to provide a cornu- 
copia of data of interest to astrophysicists and general 
relativists. 

This promising astrophysical scenario is a major moti- 
vator for much of the effort in this problem (certainly for 
the present authors). Reliable theoretical models of the 
inspiral waves will be needed both in order to maximize 
the science return from LISA data, but also to assure 
detection of these events. Ideally, the final science mea- 
surement of a detected signal will be performed with a 
"measurement template" — a model waveform computed 
accurately enough that it remains in phase with a frac- 



tional accuracy of roughly 1/p (where p is the measured 
signal-to-noise ratio) over the entire inspiral. The tech- 
niques which we describe here cannot construct wave- 
forms accurate enough for this task; as we describe in 
Sec. II Gl measurement templates require a more accurate 
and complete analysis of the small body's self-interaction 
than our techniques encompass. 

Our goal instead is to develop waveform models that 
are sufficiently reliable that they may be used to develop 
data analysis techniques for EMRI detection. EMRI 
waveforms are characterized by 14 parameters 2 . The 
number of measurement templates that would be needed 
to fully cover the 14 dimensional manifold of waveforms 
would render any search for these waves by this method 
infeasible |17|. Detection will instead be done hierarchi- 
cally using (relatively) short segments of the EMRI signal 
[l7| . Each short segment need only match coherently to a 
model for about 10 3 — 10 4 orbits. The accuracy require- 
ments on such "detection templates" are less stringent 
than those needed for measurement templates. 



B. Sketch of our calculation 

The key pieces of our formalism have been summarized 
in depth in previous papers, particularly |lj; here, we 
provide a very brief sketch largely to set the context for 
the following discussion. 

We use the Teukolsky equation 22] to calculate the in- 
fluence of a perturbation to a black hole spacetime. This 
equation describes the evolution of a complex scalar ip^ 
which is constructed from the Weyl (vacuum) curvature 
of the spacetime; it is basically a wave equation for the 
Weyl curvature, linearized around the Kerr background 3 , 
with a source. Schematically, the Teukolsky equation is 
of the form 

P 2 V>4 = T[z(t)} . (1.1) 

The general form of the operator I? 2 and the source func- 
tion T can be found in |22j ; z{t) represents the worldlinc 
of the orbiting body. Coordinate time t (which amounts 
to time as measured by distant observers) is a particu- 
larly convenient parameterization for our purposes. De- 
rived forms of Eq. (|l.l(l relevant to our analysis are given 
in Sec. IrTlAl 



2 A useful counting of these parameters is as follows: the 2 masses, 
the black hole spin S, the initial relative position rrj, the initial 
relative velocity vo, and the binary's position relative to the 
observer R. Each vector has 3 components, for a total of 14 
parameters. Including the smaller body's spin s raises the count 
to 17; fortunately, s can be neglected for our purposes llSl . 

3 Indeed, Michael Ryan has shown that the Teuk olsky equation can 
be derived from the "Penrose equation" l23l . a nonlinear wave 
equation for the Riemann cur vat ure tensor that is constructed 
from the Bianchi identity; see |24| . 
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A key thing to note at this point is that the source 
depends on the orbiting body's worldline z(t). In order 
to construct T, we use the "zeroth" order worldline, ne- 
glecting radiation. This worldline is built from geodesies 
of the background black hole spacetime; as we discuss in 
Sec. II CI we incur an important cost due to this setup. 
From this worldline, we build the source function, and 
then compute ip&. This complex scalar completely en- 
codes the GW flux to distant observers, and down the 
black hole's event horizon [25| (equivalently, the tidal in- 
teraction of the hole with the orbiting body H^|). All 
of the quantities which we use to describe GW emission 
and GW induced orbit evolution are encoded in tp^. 

We take advantage of the fact that the Teukolsky equa- 
tion is separable — ipi can be expanded into Fourier and 
spheroidal harmonics, allowing it to be computed mode- 
by-mode. Each mode is characterized by a spherical- 
harmonic-like integer index and by three integer indices 
(to, k, n) that label harmonics of the three fundamen- 
tal frequencies that describe the orbits. Thanks to the 
linear nature of the Teukolsky equation, the modes are 
independent of one another. Codes that solve for 04 in 
this manner are thus easily parallelized — many mode 
contributions can be computed separately and indepen- 
dently. Indeed, we have found that the code developed 
for this work exhibits almost perfect 1/A^ computation- 
time scaling (where N is the number of processors) 

We incur an important cost by expanding ip4 in modes: 
expanding formally requires that we understand our 
source's behavior for all time, — oo < t < oo, in order that 
the Fourier integral exist. In practice, this means that the 
orbit cannot evolve "quickly" : we require the radiation 
to backreact adiabatically, so that over a typical "orbital 
time" T or b the change in any parameter that should be 
constant is much smaller than the parameter itself (e.g., 
TorbE <C E). In principle, the adiabaticity requirement 
could be circumvented by working in the time domain 
— solving the wave equation for 04 directly rather than 
expanding in modes. Indeed, time domain codes have 
proven superior to frequency domain codes for analyzing 
problems without sources (i.e. black hole ringdown), and 
are a more natural approach to evolvin g ra diation prop- 
agating in a black hole spacetime [27], l28t 12^. l3(il ] . Al- 
though new techniques may eventually change this story 
31], at present, frequency domain codes appear better 
suited to handle the point-like sources appropriate to the 
EMRI problem. 

The next natural step is to approximate the inspiral 
and associated waveform from a sequence of geodesic or- 
bits. Beginning with some starting orbit, we compute '04. 
From it, we extract a snapshot waveform, which instan- 
taneously approximates the true waveform, as well as the 
rates of change of the orbital constants. These rates of 
change tell us how the trajectory evolves from the present 
orbit to the next orbit in the sequence. By repeating this 
process many times, we can build an adiabatic waveform 
that usefully approximates the true waveform. We have 
not yet taken this step for generic orbits, though most of 
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TABLE I: A sketch of the history of this technique. The 
magnitude of the massive black hole's spin angular momen- 
tum is given by aM. Due to spherical symmetry all geodesies 
of Schwarzschild (a — 0) black holes are planar, and can be 
considered equatorial. 

the tools needed to do so are now in hand H, IS El IS- 

Table [I] summarizes the history of this approach. To 
date, the program has been completed only for non- 
spinning black holes circular-equatorial orbits 
circular orbits j^fj, and equatorial orbits 4 0, EH- For 
the astrophysical EMRI problem, these are unrealistic 
restrictions. Since extreme mass ratio binaries are cre- 
ated through capture processes, we expect inclination 
to be randomly distributed. We also expect the eccen- 
tricity to be substantial — although radiative backre- 
action strongly reduces eccentricity, EMRI events have 
such large initial eccentricities that an estimated 50% of 
all observable EMRIs will have eccentricities e > 0.2 as 
they approach their last stable orbit [l^j. We must un- 
derstand generic EMRIs in order to realize the event rates 
predicted in Ref. |l7j . 



C. Applicability and limitations of our approach 

The "flux balancing-adiabatic evolution" approach we 
advocate here is, as emphasized above, an approxima- 
tion to the evolution of extreme mass ratio binaries. A 
more rigorous approach, upon which most workers in this 
field are focused, is based on computation of the self 
force — the small body's self interaction with its own 
gravitational perturbation. The gravitational self force 
is analogous to the Abraham-Lorentz-Dirac electromag- 
netic self force [^3, El • Its complete general relativistic 
formulation was worked out in detail by Mino, Sasaki, 
and Tanaka (3^ and by Quinn and Wald |4j|; Poisson 
41] provides a very readable summary of this formalism. 
Developing this approach to the point that one can build 
a code around it to study the evolution of generic or- 
bits of Kerr black holes remains some time in the future; 
however, efforts are intense and progress is rapid |4^ . 



4 Although Refs. 0] and l35l do not use their snapshot data to 
compute model inspirals, it is a straightforward extension to do 
so. 
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At least heuristically, the self force can be broken into 
two pieces: a "dissipative" force which carries energy 
and angular momentum away from the binary (and thus 
drives the smaller body's inspiral), and a "conservative" 
force which does not. These forces are computed relative 
to geodesies of the background spacetime — the space- 
time of the large black hole. The conservative force tells 
us that, even in the absence of radiation emission, the tra- 
jectory is modified relative to the background geodesies 
— the spacetime is deformed from that of a black hole, 
so that black hole geodesies do not precisely describe the 
motion of the small body. 

Our approach totally neglects the conservative self 
force — by construction, we can only analyze dissipa- 
tive effects. This can be seen in the schematic descrip- 
tion of our calculation near Eq. — the worldlinc 
used in the source function T is built from geodesies of 
the background spacetime, without incorporating the in- 
fluence of the conservative force. We have argued pre- 
viously that this should be adequate for scoping out 
EMRI waveforms and exploring issues related to LISA 
data analysis. Our argument was based in part on how 
the phase error resulting from our approach scales with 
the mass ratio fi/M [Eq. (1.2) of Ref. Q]. That scaling 
was in turn based on the idea that (quoting from Ref. 
@) " Dissipative terms accumulate secularly; conserva- 
tive terms do noP . In other words, dissipative aspects 
of the self force would accumulate phase effects over an 
EMRI event; we expected that conservative terms would 
oscillate, and thus not contribute as strongly over an in- 
spiral. 

To our chagrin, it has now been clearly demonstrated 
that our claim that conservative effects do not accumu- 
late secularly is wrong (although the scaling of phase er- 
rors with mass ratio which we derived from this is cor- 
rect). Using a compelling toy model to describe the influ- 
ence of the conservative self force, Pound, Poisson, and 
Nickel show that there is a component to the phase 
evolution missed by a "dissipative only" evolution. A 
simple way to understand this additional component is 
as follows: The background geodesies are characterized 
by oscillatory motion with three frequencies, fl^, fig, and 
f2 r , where fl x describes oscillations (or orbits) associated 
with the coordinate x. A conservative self force changes 
the "potentials" that determine orbital motion, and thus 
modifies these three frequencies: 

ri^ — ► n<^ + o^iff, 
n 9 -> n 9 + sn 8 

n r -> n r + sn r . (1.2) 

The shifts are smaller than their frequencies by (roughly) 
the mass ratio: 8Cl x ~ (/i/M)Q x , where \x is the mass of 
the smaller body and M is the mass of the black hole. 

Some of the most interesting strong field features are 
due to beating between these frequencies. For example, 
periastron precession (well known in the solar system due 



to planetary perihelion precession) occurs at 

fipp = — . (1 .3) 

Another effect is Lense-Thirring precession, the rotation 
of the orbital plane due to frame dragging; it occurs at 

LT = - fig . (1.4) 

The conservative self force shifts these precessions: 

<5f2pp = Sfl^ — 5Q r , 

SQur = 5^ - 6fl e . (1.5) 

It is likely that SQ^ ~ 5 fig; they are presumably exactly 
equal for Schwarzschild black holes (spherical symmetry). 
However, 5fl r is likely to differ quite a bit from the other 
two frequency shifts, importantly modifying periastron 
precession. We thus expect that eccentric orbits in par- 
ticular will be impacted by the conservative self force. 
This, indeed, is what is found by Pound, Poisson, and 
Nickel. 

Our inability to incorporate this conservative effect 
into our analysis limits its utility. We advocate our ap- 
proach primarily to begin exploring the space of EMRI 
waveforms in preparation for EISA's EMRI data analysis. 
The best understanding at the moment, after Ref. is 
that the phase errors scale like (p/M)° and (v/c)~ 3 ; con- 
sequently, in the regime v ~ c relevant to LISA the errors 
are formally of order unity. It is thus not yet clear if these 
errors are large enough to prevent adiabatic waveforms 
being useful for search templates. At any rate, though 
adiabatic waveforms may miss an important contribution 
to the phase evolution of extreme mass ratio binaries, 
they accurately represent the spectral spread that can be 
expected. Experience from circular jlll36|| and equatorial 
0, studies shows that waves from strong field orbits 
radiate significant power into high harmonics of the or- 
bital frequencies. The LISA datastream is expected to 
contain many simultaneous confused sources — certainly 
millions of white dwarf binaries (whose radiation is es- 
sentially monochromatic), perhaps ~ 10 3 simultaneous 
EMRI sources, all lying under the signals from massive 
coalescing cosmological black hole binaries. Learning to 
detangle these many signals will require models for the 
different waves which accurately describe their time and 
frequency overlap. By modeling the frequency spread of 
EMRI waves with good accuracy (though not necessarily 
the frequency evolution), adiabatic waveforms should be 
a very useful tool. 

If it turns out to be possible to separate conservative 
and dissipative self force contributions, it should be easy 
to modify our approach to include the leading conserva- 
tive effects: we would simply replace the geodesic world- 
lines presently used with worldlines that are augmented 
by these conservative effects. We would then use this 
augmented worldline to construct the Teukolsky source 
function. Doing so would allow us to calculate an inspiral 
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that accurately accounted for the leading dissipative and 
conservative effects 5 . 

It should be strongly emphasized that this is a rather 
big "if": Because the gravitational self force is a gauge 
dependent quantity (though its impact on binary orbits 
must be gauge independent), it is far from clear that one 
can cleanly separate the "conservative" from the "dis- 
sipative" influence. It may be that the only way to dis- 
criminate the two influences is to run the equations of mo- 
tion which the self force implies through a self-consistent 
wave-generation calculation |43| . in which case there is 
no need to separate these influences in the first place. 



D. Organization of this paper 

We begin this paper by summarizing, in Sec. ^ rele- 
vant properties of bound Kerr black hole orbits. These 
orbits determine, in turn, the properites of the Teukol- 
sky equation's source term. Our goal is to present enough 
detail to make it clear how we build this term. Of par- 
ticular importance is a frequency domain des crip tion of 
functions that are built from these orbits [44l l45| . 

We next (Sec. IIII|) summarize the relevant details of 
the black hole perturbation formalism that we use. We 
go into some detail in this section. This is in part to 
make the analysis self-contained, but also to correct some 
small errors that have appeared in previous papers, no- 
tably Ref. (hopefully without introducing too many 
new errors). We first review in Sec. llll Al how Teukolsky 
equation solutions are built. This is a somewhat sub- 
tle numerical problem — due to the long-rangedness of 
the separated Teukolsky equation's radial potential, the 
amplitudes of ingoing (oc e~ %ulr ) and outgoing (oc e %U3T ) 
radiative solutions grow at different rates. [The quantity 
r* is the Kerr "tortoise coordinate" ; cf. Eq. i|3.16|l below.] 
One component of the solution can easily swamp the 
other, causing numerical resolution to be lost. Sasaki and 
Nakamura (46J found an equation with short-ranged po- 
tential whose solutions are related to those of the Teukol- 
sky equation by a simple differential transformation; the 
combined Sasaki-Nakamura- Teukolsky formalism makes 
for a very well-behaved numerical problem. We next 
(Sec. IIII Bfl discuss explicitly how the source function is 
built, taking advantage of the results that we presented 
in Sec.[n] and then discuss some pratical issues related to 
the numerics fSec. IIII^ . We wrap up this section with 
a brief discussion of how the gravitational waveforms and 
fluxes of "conserved" quantities are extracted from these 
solutions. 

Section IIVI discusses several practical issues related to 
numerics. We first describe the algorithms used to nu- 
merically represent the geodesies, a key element for the 



5 One might label a self force that allows this separation "compas- 
sionately conservative" . 



Teukolsky source term, before describing issues related to 
computing each of the modes from which we build ip^. We 
then discuss in great depth the algorithms we have devel- 
oped to truncate our modal expansion. Strictly speaking, 
the number of modes that should be used to build ipi is 
infinite; picking a finite value to truncate this expansion 
that is "large enough" is somewhat subtle. We describe a 
scheme in Sec lIVCl that has proven to be robust enough 
to work well for our present purposes. There is clearly 
room for improvement, however; we describe some ways 
in which this procedure could be made better. Finally, 
we conclude in Sec. lIVDl with detailed discussion of vali- 
dation tests that we made against previous results in the 
appropriate limits. We find that the code agrees per- 
fectly with results from Ref. 1] when the eccentricity is 
zero 6 , and agrees quite well with the "equatorial" code 
of Glampedakis and Kennefick in Q in the limit of zero 
inclination. Except for some (very small) down-horizon 
modes that do not significantly impact the overall flux, 
we typically find agreement at the level of 10 -3 or bet- 
ter with Glampedakis and Kennefick. This is the level of 
accuracy that they cite for their code. 

Our main results are presented in Sec. E] The modal 
decomposition allows us to break the waveform into 
"voices" : 

H = h+-ih x =J2 H kne- iUkn(t - r ' ) , (1.6) 

fcn 

where LOkn = kflg -\-nVL r . [A sum over / and m, including 
oscillations with frequency mfi^, is hidden in the defini- 
tion of Hkn', see Eq. (|5.3|l below.] The "polar voice" is 
composed of the terms with n = and k ^ 0; The "radial 
voice" is composed of the terms with k = and n ^ 0; 
The "azimuthal voice" is the term with k = n = 0, and 
the "mixed voice" is composed of the terms with k ^ 
and n ^ 0. A similar voice-by-voice labeling can be ap- 
plied to the fluxes. We find that the importance of the 
various voices is fairly simply controlled by the eccen- 
tricity and the inclination angle: Polar voices become 
progressively more important as the purely polar orbit 
(inclination 90°) is approached from either direction; ra- 
dial voices rapidly become important as eccentricity, e, 
approaches unity, typically dominating the waves' spec- 
trum by e ~ 0.3 — 0.5. We speculate that this multivoice 
character may facilitate approximations in the design of 
GW detection schemes, making it possible to detect the 
most important voices of signals, rather than needing to 
detect the (rather ornate) chorus of all voices together. 

Section IVT1 presents our final conclusions, and outlines 
future work to which we plan to apply this formalism. 



6 The "generic" code is a direct descendent of the "circular" code 
used in Ref. so it is perhaps not too surprising that there 
is agreement in the circular limit. However, many details of the 
generic code's inner workings are significantly different from the 
circular incarnation, so this agreement is far from trivial. 
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Chief among these future tasks will be augmenting this 
approach with an adiabatic scheme to evolve the Carter 
constant 0, 0- an( A then using the complete flux data 
to build model adiabatic inspirals and their associated 
waveforms (as discussed in Sec. IIB|I . We are optimistic 
that this can be completed relatively quickly since there 
appear to be no major hurdles or issues of principle that 
must be overcome first. Producing adiabatic waveforms 
for initial data analysis algorithm development will then 
"just" be a matter of finding sufficient CPU power. 



II. BOUND BLACK HOLE ORBITS 

In this section we review the geodesic motion of a non- 
spinning test mass on a bound orbit of a Kerr black hole. 
Kerr orbits are not a new subject of investigation |3.l47|. 
but interest has been renewed recently because of their 
relevance to EMRI GWs H0IH. Thc 

main new de- 
velopment has been Mino's exploitation of the fact that 
these orbits are fundamentally periodic entities 0- The 
utility of exploiting this property is discussed in detail in 
Ref. |44j , and details of its application to EMRIs can be 
found in Sec. 3 of Ref. Q. Here we will discuss the rela- 
tion between an orbit's geometry and its natural frequen- 
cies. Using a special choice of initial conditions (the fidu- 
cial geodesies from Ref. Q) we also derive a new system 
of equations for efficient numerical evaluation of these 
orbits. 

Kerr black holes are characterized by two parameters: 
the mass M of the black hole, and the magnitude aM of 
its spin angular momentum, with < a < M. Through- 
out this paper, we will use Boyer-Lindquist ^3 coordi- 
nates (t,r,9,4>), with units chosen so that G = c = I. 
The line element for the Kerr geometry is then given by 



3] 



a ">Kcrr 



dt 2 + — dr 2 + S dO 2 
A 



2 i 2 

r + a 



2Ma 2 r 



d<t> 2 



AMar 



where 



r 2 + a 2 cos 2 I 



■ sin 9 dt d<f>, 



2Mr + a 2 



(2-1) 



(2.2) 



Bound black hole orbits admit four constants of the 
motion which allow us to rewrite the geodesic equations 
as a system of first order differential equations. Three of 
these constants are fairly straightforward — if the test 
particle has 4-momentum p, these constants are 



described by a system of first order equations . As dis- 
cussed in the Introduction, it is often useful to think of 
the Carter constant as representing the geodesic's non- 
axial angular momentum (a correspondance which is ex- 
act for non-rotating black holes). 

Carter's first order equations can be written in the fol- 
lowing form 



2 

= V e (0), 



dt 
dX 

d(t> 
dX 



V t (r,6), 



(2.4) 



The functions Vt, V r , Vg, and V$ are shown in Appendix 
1X1 The Mino time parameter A is related to the test 
mass' proper time r by 



dr_ 

dX 



S. 



(2.5) 



Mino time decouples the radial and polar equations of 
motion so that V r = V r (r) and Vg = Vg(6). This prop- 
erty appears to have been recognized by Carter; however, 
Mino appears to be the first to use A to fully exploit the 
periodic nature of Kerr orbits Since the first order 
equations for r and 9 are purely quadratic in their deriva- 
tives, r(A) and (9(A) are periodic functions for bound or- 
bits. The functions dt/dX and d<j>/d\ are then biperiodic 
functions in the sense defined by Ref. ^4[. As a result, af- 
ter subtracting a term proportional to A, the coordinates 
t{X) and 4>(X) can be represented with a two dimensional 
Fourier series. These properties, discussed in more detail 
below, are extremely powerful. 

Solutions of the geodesic equations 1)2. 4f) are uniquely 
determined if we specify E, L z , Q and thc initial position 
z(X = 0) (or some other equivalent set of constants). We 
now focus on the orbit's geometry or shape, which is 
determined by E, L z , and Q. This fact can be roughly 
understood in that the orbit must be bounded by two 
radii r m - m < r < r max , and [because the geometry (|2.f I) 
is symmetric across the equatorial plane] one polar angle 
#min < 9 < (n — 9 min ). The orbit is then confined to a 
toroidal region sketched in Fig.^ The boundaries of the 
orbital torus could have equivalently been determined by 
an eccentricity e, a semilatus rectum p, and an inclination 
9 U1C defined by 



pM 



pM 



(sgni 2 



n 2 ■ 
(2.6) 

Here we have included a factor of sgnL z so as to make 
#i nc most closely resemble another common definition for 
the orbital inclination angle i 0: 



p-p = -(i 



d t -p = -E, d^-p = L zl 



(2.3) 



where \i is its rest mass, E is the orbit's energy, and L z is 
its axial angular momentum. Carter discovered a fourth 
constant Q which allows the motion to be completely 



cos L 



(2.7) 



A nice property of the angles 9- mc and i (not shared by the 
angle 9 m { n ) is that they automatically encode a notion of 
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FIG. 1: The orbital torus and the evolution of r and cos# for a generic geodesic orbit. The magnitude of the black hole's spin 
is a — 0.998M. The orbit shown here has eccentricity e = 1/3, inclination Sine = 30°, and semilatus rectum p = 7. We start 
the orbit at (r, cos 6) = (r m i n , cos # m in), and end it around (r, cos (9) = (10.5M, — 0.1), after several complete radial and polar 
oscillations. The orbit is not closed: over time, it would eventually fill the orbital torus. 



prograde and retrograde — prograde orbits (0 motion 
parallel to the hole's rotation) have 9 mc ,i < 90°; ret- 
rograde orbits (<f> motion antiparallel to the hole) have 
#inc,t > 90°. We have found that in general l ps 9 inc . 
Even in the strong field the two quantities typically dif- 
fer by less than 10%. For the orbits which we discuss 
in detail here (Sec. i and 9- lric differ by no more than 
0.5%. 

Explicit algebraic relationships between the geometric 
orbital parameters (e,p, ft nc ) and the physical constants 
(E, L z , Q) were first computed by Schmidt 0, and are 
shown below in Appendix [S] We find that when explor- 
ing the orbital parameters space, it is best to first think in 
terms of (e,p, 9i nc ), and then to convert these (following 
Appendix lA"|l to (E, L z , Q) which are used in all further 
calculations. 

Using a special choice of initial conditions, we now de- 
rive a set of equations for efficient numerical evaluation 
of bound orbits. Since the solutions to the radial and 
polar geodesic equations (|2.4[l are periodic in Mino time, 
each can be expressed as a Fourier series, 

oo oo 

9(X) = J2 0ke- lkTe \ r(A)= ]T r n e~ m ^\ 

k= — oo n= — oo 

(2.8) 

where 9k and r n are constants, and where T rj g are the 
orbital frequencies in Mino time. In Appendix lAl we list 
explicit expressions for these frequencies as functions of 
the orbital parameters described in the previous section. 
It is convenient to write these series in terms of angle- 
variables, 



The result is 

oo oo 

6(w e )= ke~ tkwe , r(w r )= £ r n e~ m ^ . 

k=-oo n= — oo 

(2.10) 

We now specialize to a specific choice of initial condi- 
tions for the radial and polar motion. We require that 
the orbit begins at a turning point of both the r and 9 
coordinates, 

r(w r = 0) = r min , 9(w e - 0) - 9 min . (2.11) 

This choice of initial conditions will result in a simplified 
numerical evaluation of the geodesies. For example, the 
coordinates r and 9 are now even functions of w r- g, 

r(-w r ) = r(w r ), 9(~wg) = r(w e ), (2.12) 

and the Fourier series 1)2. 8f) become cosine series, 

oo 

9(w e ) = 9 + 2Y / 9 k cos(kw e ), (2.13) 

k=l 
oo 

r(w r ) = r + 2 r n cos(nw r ). (2-14) 

ri=l 

In this paper we will not evaluate these series explicitly 7 , 
though we will use the fact that they are even functions. 

We now derive a similarly simplified series expansion 
for t and 4>. Both the t and cf> coordinates will be treated 



w r = T r A, we — TgA. 



(2.9) 



7 Current investigations suggest that doing so may lead to greater 
computational efficiency. 
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in a similar way, so to save space we define 
x = t,<j>, x = V u V ( f n 



(2.15) 



such that an overdot represents a derivative with re- 
spect to Mino time A. Followi ng t he analysis preceeding 
Eqs. (3.18) and (3.19) in Ref. |44|. we write the deriva- 
tives x in the form 



i(A) 



OO 

E 

k=l 



(4e 



-ikTeX 



+ C.C.) 



n=l 



— inT r A 



+ C.C.), 



2tt 



/ dwg x e [9(wg)]e tkws 



1 



2n 



— / dw r x' \r 
2tt 



{w r )]t 



(2.16a) 

(2.16b) 
(2.16c) 



Here "c.c." means the complex conjugate of the preceed- 
ing term. The leading constants are defined by 



±00 = x + Xq . 



(2.17) 



We have made use of the fact that the derivatives of 
both t and <f> separate into a sum of two functions, each 
depending only on one of the coordinates 8 , 



x = x e {9)+x r {r). 



(2.18) 



Following the notation of Ref. [44j , we use the following 
symbols for x m in the cases x = t and x — <p: 



r = (v t ) 



00. 



t = m 



00- 



(2.19) 



The constant T is the analog of the Lorentz factor of spe- 
cial relativity. For example, in the case of an orbit which 
is both circular and equatorial, we have T = dt/dX = 
(dt/dr)T, (E is constant for circular-equatorial orbits). T 
also relates the Mino time frequencies g >r to coordi- 
nate time frequencies f2 r , fig, and Qd>: 



O - T * 



u e - 



T 

O = — 

"r - p • 



(2.20) 



Schmidt ^3 provides an elegant derivation of closed form 
expressions for the orbital frequencies. His results were 
converted into the Mino time frequencies in Ref. 01 ( see 
Appendix El for details). 

As a consequence of their initial values (|2.11|l , the func- 
tions 9(wg) and r(w r ) are even. The functions x e (wg) 
and x r (w r ) are then also even, and the series (|2.16|) above 



It's interesting to note that Teukolsky has shown that the master 
equation 13.21 separates in any coordinates where this property 
I2TT51 is preserved 



simplifies to 



x(X) = xqq + 2 x 6 , cos(fcTgA) + 2 x r n cos(nT r A). 



k=l n=l 
-r'[ = — I dwg X 9 [6(wg)] COs(kwg), 

dw r x r [r(w r )] cos(nw r ). 

If we now assume initial values for t and 
t(X = 0) = , </>(A = 0) = , 





1 


x k — 






7T 




1 


• T 






7T 



(2.21a) 
(2.21b) 

(2.21c) 



(2.22) 



so that we are now discussing the fiducial geodesies from 
Ref. 0, a similarly simplified expression for the coordi- 
nates x can be found by integrating the series l|2.21af) for 
their derivatives x: 



x{\) = x 00 X + Ax e (X) + Ax r (X) 



Ax°(X) = 2^ A4 sin(fcT e A), 
fe=i 

OO 

Ax r (X) = Ax n sin(nT r A), 



(2.23a) 
(2.23b) 

(2.23c) 



Axl 



Axl 



n=l 
2 



2 



mrT T 



dwg x (wg) cos(kwg), (2.23d) 
dw r x r (w r ) cos(nw r ) . (2.23e) 



Note that in the cases of orbits with special geometries, 
the coordinates x simplify accordingly: 

Ax r (X) = = Ax r n (for circular orbits) , (2.24a) 
Ax e (X) = = Ax d k ( for equatorial orbits) . (2.24b) 



Evaluating the integrands in Eqs. I|2.23dj) and 12.23c|) 
as functions of r(w r ) and 9(wg), would first require a 
direct integration of the radial and polar geodesic equa- 
tions. Unfortunately this is somewhat difficult because 
the derivatives V r and Vg vanish (by definition) at the or- 
bital turning points. This means, for example, that the 
integrand in the direct expression for r(A), 



A = 



p W dr' 



(2.25) 



contains singularities which, although integrable, com- 
plicate the numerics. This is a well known problem (see 
the Appendix of Ref. 0] for more details) which can be 
avoided if we change the integration variables to more 
well behaved coordinates x(^) an d ip(r), where 



cosx 



cos 9 
cos 9 mu 



pM 



1 + e cos ip 



(2.26) 
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Recall that p is the semilatus rectum of the orbit, and 
e is the orbital eccentricity. So in practice we replace 
Eas. l|2.23d|) and (|2.23e|) with 



dwg 



Axl = — x , 

k-Kie J ax 



Axl 



mrT r 



d X —x a (we)cos{kw e ), (2.27) 

a>X 

dtp ^Yj^x r (w r ) cos(nuv). (2.28) 



dip 



In order to use this result, we of course need explicit ex- 
pressions for wg{x), W r {ip)i and their derivatives. These 
functions were first derived in Ref. ^ij, and the results 
are shown below in Appendix 1X1 It turns out that w${x) 
and dwg/dx can be written in terms of elliptic integrals. 
The functions w r (ip) and dw r /dip are slightly more in- 
volved. 

We have found that this technique is an extremely effi- 
cient way to evaluate bound orbits and functions of them 
(see Ref. |44| for a simple example) . The two expansions 
given by Eqs. I|2.23b|l and (|2.23c[l converge very rapidly. 
Typically only about 50 terms are needed in order to 
obtain a fractional accuracy of 10~ 12 . 



III. PERTURBING A BLACK HOLE WITH AN 
ORBITING TEST MASS 

In this section we will discuss how a non-spinning test 
mass /x on a bound geodesic perturbs the Kerr geome- 
try (|2.ip . The radiative information describing the lin- 
ear order perturbation can be extracted from the Weyl 
curvature scalar 49] 



■04 = —C a p 1 &n a mPn 1 m & 



(3.1) 



Here, overbar denotes complex conjugation, C a b c d is the 
Weyl curvature tensor (the Riemann tensor in vacuum), 
and the vectors are elements of a Newman-Penrose null- 
basis [S^OOl {l a , rn a , fh a ,n a }. For distant observers, the 
curvature scalar -04 is simply related to the metric pertur- 
bation. For these same observers, and also for observers 
at the event horizon, 04 is simply related to the fluxes 
of energy and angular momentum (see Sec. IIII Dl for de- 
tails). Since we are only analyzing perturbations to first 
order in the mass ratio, 04 always represents the leading 
0(fi/M) contribution to the curvature perturbation. 



A. Teukolsky-Sasaki-Nakamura formalism 

Tcukolsky showed that 04 is a solution to an equation 
of the form [H 
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-7. 



(3.2) 



where Tis the source term (described below in Sec. IIII Bl . 



<P = P 04, 



-(r — iacost 



(3.3) 



and where /7t^ r (r) and Vt<pg(9) are second order differen- 
tial operators containing derivatives with respect to the 
variables shown as subscripts [see Eqs. IjBlall and ljBlb|) ]. 
Equation l|3.2l) is known as the master equation. In the 
next two sections, we will summarize Teukolsky's tech- 
nique for solving the master equation by separation of 
variables (2^, converting Eq. I|3.2|l into a pair of ordi- 
nary differential equations. This technique is referred 
to as solving the master equation "in the frequency do- 
main". As briefly discussed in the Introduction, one 
could instead solve the master equation numerically as 
a partial differential equation "in the time domain" (see 
Refs. H^E^I and references therein). Such an approach 
can be advantageous for problems where there is no 
source, or where the source is not pointlike. To date 
at least, time domain treatments are accurate to about 
10% [53| for EMRI sources (at least for quantities like 
orbit averaged fluxes of energy and angular momentum) . 

In the source free case T = 0, the master equation l|3.2|) 
is satisfied by functions of the form 

<p m (u,C) = R m (r;uj,C)S m (9;uj,C)e~ M "^, (3.4) 

where C is some constant, m is an integer (since ip must 
be periodic in 0), and oj is any real number. The func- 
tions R and S are solutions to ordinary differential equa- 
tions of the form 



U r (r,m,uj) - C R m (r,u),C) = 0, 
V e (e,m,u) + C] S m (9,uj,C) = 0. 



(3.5a) 
(3.5b) 



Requiring that the solutions to the angular equation 
(|3.5b|) (known as spin-weighted spheroidal harmonics 
with spin weight —2) be regular results in a discrete spec- 
trum of eigenvalues C = Ci m (u), 



S m [Q,uj,Ci m (uj)] = Si m {0,uj) , 



(3.6) 



where / > max(|m|,2) is an integer. The functions 
Si m {6,uj) are uniquely defined only after we specify 
boundary conditions, or equivalently after we chose a 
normalization convention. We use the convention 9 



/ d9 {S lm (d,uj)} 2 S m9= ±- 



(3.7) 



The choice of normalization is arbitrary in the sense that 
it does not change any physical predictions. However, in 
order to compare equations from other papers in which a 
different convention was used, we must state our choice 
explicitly. The functions Si m (9,u) will be used as a basis 
to express functions f(t, 9, 0) as follows: 



1=1 



a— — l 



duo {lmuu\f)S lm {9 : u)e- wt+ml 't' , 

O 

(3.8) 



9 Although not stated there, this is the normalization convention 
used in Ref. [lj. 
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where 



(Zmu;|/) = i Jdt Jdn f(t,6,cb)S, 



(3.9) 

so that we have {V m! uj 1 \lmuj) = 5u>5 mm i8(u} — lj'). 

We now return to the case where the source T in the 
master equation Ij3.2|l is nonvanishing. Assuming a solu- 
tion of the form 



OO I r, 

( O ^ 1 J 



(3.10) 



1=2 m=-l ' 

where 

<p lm (w) = Ri m (r,u)Si m (0,u)e- Mm *, (3.11) 
implies that the radial function i?; m (r, u>) must satisfy 
d ( 1 d 



dr \ A dr 



Vi m {r,uj) R lm (r,uj) = -7 lm (r,u), 
where T; m (r, u>) = (lmw\7), and the potential is given by 



K 2 + 4i(r — M)K 



Here K = (r 2 + a 2 )uj — ma, and 

Az m (w) = -Ci m (u) - 2amuj. 



+ 8iu>r + \i m (u). 

(3.13) 



(3.14) 



The radial equation l|3.5a|l is just the source free (homo- 
geneous) version of the previous radial equation i|3.12|l . 
The two independent solutions to the homogeneous ra- 
dial equation will be used to construct the solution to 
the inhomogeneous equation. We use the so-called "in- 
up basis" {Rjl n (r, oj), Rf^irAv)}, which is defined by the 
following limiting behavior p|: 



Rhn(r -» r+,u) = B; h „° i l0 (w)A 2 e 



2-iPr* 



(3.15a) 



Rf m {r -> oo,w) = BZ\u)r 3 e^' + -i^e^', 

(3.15b) 

noo /„ , „ , A r-iout / \ iPr* , nin / \ A 2 —iPr* 

(3.15c) 
(3.15d) 



R^(r - oo,w) = D^ye^r 



Here P = uj — ma/(2Mr+) ) and the "tortoise coordinate" 
r* satisfies dr* /dr = (r 2 + a 2 )/ A: 

2Mr + r — r + 2A/r_ r — r_ 
r (r) = r H In — r- In ■ 



r+ - r- 2M r+ - r_ 2M ' 

(3-16) 

where r± = M ± v 7 M 2 — a 2 are the roots of A (r = r+ 
is the location of the event horizon). 

The numerical calculation of the homogeneous solu- 
tions, say iifL(r), would in principle work as follows. 



First set Bf^uj) = 1 (a normalization convention). Then 
note that 



(3.17) 



(3.18) 



Next starting near the horizon at r = r + , integrate for- 
ward. When we reach r s» oo, we can read off i3 out /i? llolc 
and l/£> llolc . (When clarity permits, we will often drop 
the somewhat cumbersome subscripts I, m and the depen- 
dence on r, oj.) Since the first term in Eq. (|3.18|l grows 
rapidly with r, the two terms cannot be extracted with 
equal accuracy. We work around this by instead solving 
the homogeneous Sasaki-Nakamura equation [4(| (which 
was designed to avoid this problem) and converting the 
result to R H o °] see Ref. [l| for details. Recently Fujita 
and Tagoshi have worked out a sophisticated numerical 
scheme for computing the functions R H '°° numerically 
with great accuracy (with fractional accuracies ~ 10~ 14 
in the corresponding energy fluxes for orbits which are 
both circular and equatorial) [54|. 

The general solution to the radial equation l|3.12(l , cor- 
responding to the retarded solution of the master equa- 
tion (|3.2|l . is 

Ri m (r,uj) = Z^(r,u})Rf^(r,oj) + Z^(r,w)R^ n (r,u), 

(3.19) 

where the functions Z are radial integrals over the source 
term: 



Al m (oj) 



A' 2 



7i m (r',u), 



(3.20a) 



Z%(r,uj) = / dr' * lm Z' h l m {r',u) 



Al m (v) J r 



The function Ai m {^>) is given by 1 



Ai m (oj) = 2iojBZ(oj)D^(oj). 



(3.20b) 



(3.21) 



This construction of Z K '°° and A follows from the theory 
of Green's functions • 

In this paper, we only evaluate the perturbed field ip± 
[and therefore the functions Z (|3.2U[) ] outside the orbital 



torus 



r < r n 



and r > 



Evaluating ^4 at arbi- 



trary locations within the orbital torus r m j n < r < r max 
would require a more complicated numerical apparatus 
than the one developed here. 



The value of A reported in Ref. [j| has the wrong sign. 



11 



B. Bound test mass sources 



where 11 157 



Beginning from Teukolsky's original expression for the 
source term T |23 , Breuer [56j has com put ed the explicit 
form of the projected source (see Ref. [57j or Sec. IV C 
of Ref. 0): 



7 lm (r,uj) = J dQdt % n {t^e^,u)S lm {9 1 u;)e wt - im ^ 

(3.22) 

where 



2, 



V2^p^L^[p- i p 2 J + {p- 2 p- 2 ^ 1 T nrfl 
-2p 3 L_ 1 [p- 4 L (p- 2 p- 1 T„ n )] 



-V2A 2 p 3 J+{p- 4 p 2 A- 1 L- 



A 2 pV+[p- 4 J+(p-^7k rft )] 



(3.23) 



This function is given in terms of the tetrad components 
of the orbiting particle's energy-momentum tensor, 



(3.24) 



where T a0 is the energy-momentum tensor of the parti- 
cle, and a and b are either norm. We define the deriva- 
tive operators 



J+ = d r + iK/A 
L« = da + esc 6 



(3.25) 

aw sin # + scot#. (3.26) 



In Kerr spacetime, a point particle of mass p has an 
energy-momentum tensor given by 



T a f>(t,r,e,<t>) = 



pu u 



6[r-r(t)]6[6-e(t)}8[<f>-<f>(t)l 



S sin oft 

(3.27) 

where u a is the particle's 4- velocity. The quantities r(t), 
9(t), and <fi(t) appearing in Eq. i(3~2Tjl are the coordinates 
of the geodesic at coordinate time t, and should not be 
confused with r, 9, and <f> in the definition of the source 
term (jS23>- 

We use the Newman-Penrose basis given by Kinnersly 
[see Eqs. (jB2j| and JB3j|]. This means the tetrad 
components of T a ° are given by 



C 

Cfh 



dX p 
d\ pp 2 



dr 



E (r 2 + a 2 ) -aL z + — 



•IP L * V - ■ d0 

i ah ^ — sm 



sin 2 9 J "' ' dX 



dt 2 

dX pp 
~d~t 2V2E 

i I — I sm 9 H 

aA 



dr 
dX 



(3.29a) 



(3.29b) 



(3.29c) 



These functions (|3.29|) depend explicitly on dr/dX and 
d9/dX, rather than on their squares V r fi. Extra book- 
keeping is thus needed to keep track of the signs of dr/dX 
and d9 / dX when evaluating them numerically. 

Following the details shown in (for example) Ref. pj, 
we can now write the projected source term H3.22(l in the 
following form 



Tz m (r,w) 



dt A 2 {[A nn0 + A nr - n0 + AnmoMf - r(t)] 

+d r [(A i)S(r - r{t})} 

+S r 2 [^ 2 <5(r-r{<})]} e ^-W . (3.30) 



The quantities A aoc are shown explicitly in Appendix [H] 
Substituting Eq. (|3.30jl into Eq. H3.20jl and eliminating 
the radial delta functions gives 

Zf m (r,u>) = -~ I dt e iut - im ^Q*[r,r(t)] 



[{A nn o + A nm Q + A mm o) 

(A A \ d A ^ 

dr ar* 



(3.31) 



where * = H, oo. We define the step functions 

e°°(xi,sa) = Q li (x 2 ,x 1 ) =Q(x 2 -xi), (3.32) 

in terms of the Heaviside step function 0(x). All other 
functions of r and 9 under the integral in Eq. (|3.31|) are 
to be evaluated at r(t) and 9{t) respectively. Because 
r only appears inside the step functions, the quantities 
Z H, °° are independent of r for all r > r max and r < r m ; n . 

To clean up this expression, we now absorb most of the 
integrand into a single function: 



dt e u 



(3.33) 



T„ 



Cab 



5[r-r{t)]S[0-0(t)]d[4>-<l>(t)], (3.28) 



11 In Ref. l5Sl . the dd/dr terms are missing from the expressions for 
Cnm and Cmfh- Also, in Ref. dd/d\ is erroneously replaced 



with (de/dX) 2 
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Following the arguments in Sec. V B of Ref. j44jj we 
exploit the harmonic structure of the geodesies (Sec. |HJ) 
to simplify Eq. (|3.33() . In the end, we will arrive at a 
general expression for ipi 1)3. 1U[) as a discrete sum over 
frequencies, rather than an integral. First we insert 



Numerical considerations 



t = TX + At, 



T A + A(/>, 



(3.34) 



into Eq. (|3.33|) . and we change the integration variable 
from t to A. The result is 



2L(r,w)= d\e^ T ~ mr ^Jt m {\r,Lo), (3.35) 



where 



dt 



J? m (\, r, oj) = —e^-^If^X, r, u) . (3.36) 

The function J* m (X, r, u>) is biperiodic, so we can write it 

as 441 



i{kT e +nT r )X 



(3.37) 



kn 



the constants Ji rakn are given by 



T* 

Imkn 



l 



27T p2lT 

dw e / dw r e ^ kw o+nw r ) 



(27T) 2 Jo 



(3.38) 



By inserting the expansion (|3.37(l into Eq. H3.35JI . we find 
that the integral 1|3.35[) is just a sum of delta functions: 

ZL(r, W) = Z Lkn(r)S(" ~ ^mkn) ■ (3.39) 
k n 

We have used the coordinate time frequencies (|2.2Q(I to 
define 

u m kn — rnQ^ + kQg + nfl r . (3.40) 
The expansion coefficients for Z* m {r,uj) are given by 

X Jf m ( w r,W0,r,UJmkn), (3-41) 

By substituting the expanded form H3.39|) of Z into the 
general expression (|3.1U|) for ^4, we eliminate the fre- 
quency integral to obtain 



i/ji = p 4 ^ Rimkn{r)Si mkn {9)e 



(3.42) 



Imkr, 



The radial and angular functions are now discrete func- 
tions of frequency, inheriting the angular harmonic index 
k and the radial harmonic index n: 

SlmkniQ) =Sl m (d,CJ m kn), (3.43) 

=Rlm(r,Umkn), (3.44) 

Rimkn(r) =Zj i nkn (r)R'j% lkn (r) + Z^ lkn {r)R^ mkn {r). 

(3.45) 



Because the functions r(w r ) and 0(wg) are periodic 
and even [cf . Eq. l|2.12Jl ] , the range of the integrals (|3.41(l 
can be reduced by a factor of two. Writing the integrand 
in (|3.41() as z(w r ,wg), we can expand the integrals into 
the form 



PIT Pit p27Z 

Z = I dwg / dw r + / dwg I dw r 

JO Jo Jit 

2tt c-ix p2tz p2tt 

+ / dwg / dw r + / dwg / dw r 

JO Jtt Jtt 



z(w r , w e ). 

(3.46) 



Since the integrand has a period of 2ir in both variables, 
we can shift all of the {ir, 2ir} branches to {— w, 0}: 



/•7T />7T />0 

Z = I dwg / dw r + / dwg / dw r 

JO JO J-TT 

pir pO pO 

+ I dwg / dw r + / dwg / dw r 

JO J -TV 



z(w r ,we). 

(3.47) 



We can now reflect all of the {— 7r, 0} branches across the 
origin to get 



/>7T />7T 

Z= / dwg / dw r z(D r W r , DgWg) 

J ° J ° D r =±D e =± 



(3.48) 



Direct evaluation of Eq. (|3.48|l is complicated by the 
fact that it is difficult to evaluate r(w r ) and 9(wg) [see 
discussion preceeding Eq. (|2.27|l ] . We avoid this problem 
by changing integration variables from w r and wg to ip 
and x [see Eq. ifT^ ]. 

Eccentric orbits introduce another numerical nuance. 
The integrand in Eq. (|3.48|l is moderately divergent if the 
orbit covers a large range in r (i.e., if the orbit is highly 
eccentric). This problem can be traced to the scaling 
behavior 1)3.15(1 of solutions to the homogeneous radial 
equation. In particular, the typical size of the integrand 
at r m i n can be different from typical values at r max by 
orders of magnitude. This causes an artificial focusing of 
numerical accuracy. Our choice of boundary conditions 
(|2.11() is such that the turning points r mm and r max occur 
at ip — and ip — n respectively. To avoid this problem, 
we introduce a change of integration variables from ip to 
£ where dip/dQ = e - ^'* , with * = H, oo, and 



CWO = ^ (e^ l) 



(3.49) 



Here '&* is a constant (which can be positive or negative) . 
Before integrating, we tune the value of so that the 
typical magnitude of the integrand at the outer turning 
point is comparable to that at the inner turning point. 
In general, the value of will decrease if either I or 
eccentricity e are increased. 
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Applying these various changes of variable and mas- 
saging of the integration ranges to Eq. I|3.4ip . our final 
expression for the Zf mkn (r) coefficients is given by 



Zlmkn ( r ) 



1 

"2nT 



D r =± D„=± 

x e - ^'* I* m (w r ,w g ,D r ,D e ,r,uj mkn ) 
x exp [iDg (kwe + Lu m k n At e - mA(f> )] 

x exp [iD r (nw r + w m fc n At r — mAcjf)] , 

(3.50) 

where I* m is defined by Eq. I|3.33[) . The functions r(t), 
9(t) and their derivatives are evaluated as 



r(i) -> r(C), 



r(t)] 



D r 



dr 
dX 
d6 
d\ 



HO] 
[0(x)} 



(3.51) 



Lastly, note that applying the simultaneous transfor- 
mations m — > —to, a; — > — a; to the radial (|3.12|) and the 
angular equations l|3.5bjl . one can derive a simple sym- 
metry rule for expansion coefficients 0, 



zr, 



f(-m)(-fc)(-n)( r ) — ( — 1) Zlrnkn( r ) 



(3.52) 



This symmetry is used to reduce the computation time 
for various summations (see Sec. II V (J)) . 



D. Waveforms and fluxes 

In this section, we use the expanded form of the cur- 
vature perturbation (|3.42|) to derive expressions for the 
leading order metric perturbation at infinity, and for the 
radiated fluxes of energy and angular momentum to in- 
finity and down the black hole's event horizon. Since we 
will be interested in quantities at both infinity and the 
horizon, it will be useful to first define the limiting values 
of the coefficients Z* mkn (r) as follows: 



7°° — R h °l c 7°° ( r ^ r ■ \ 
^Imkn ~ ^ 'irnkn^ 'ImknV ' 'mm)- 



(3.53a) 
(3.53b) 



Recall that the functions Z* mkn (r) are independent of r 
for all r > r max and r < r min . 

At infinity, the leading order curvature and metric per- 
turbations are related by 

V> 4 (r^oo) = ^(/ l+ -i/ lx ). (3.54) 

Here h + and h x are the two independent components of 
the metric perturbation, defined by 

h aP {r -» oo) = h+e+g + h x e* + 0{l/r 2 ) , (3.55) 



where e\ B and e*^ are polarization tensors [6C 



If we 



now substitute the expanded form of the curvature per- 
turbation H3.42(l into the left hand side of Eq. H3.54jl . 
evaluate the limit, and integrate twice in coordinate time 
(setting the arbitrary integration constants to zero), we 
find 



9 7 H 
£ \ f^rofcrt 

r ^ J 1 , 

Imkn mkn 



Slmkn{d)e 



-ioJmkn {t—r)+im4 



(3.56) 



We have used the definition l|3.53a|) for the limiting value 
of^(r). 

It is important to note that the expression for Zj^^ 
derived above [e.g. Eq. I|3.50|l ] assumed geodesies with 
a specific choice of initial position: r(0) = r m ; n , 8(0) = 
Omin, t(0) — 0, and </>(0) = 0. These are the fiducial 
geodesies from Ref. Q. Making a different choice of ini- 
tial position will result in an overall phase change for 
^imkn • The details of this phase can be found in Ref. • 
As we will see below, the expression for the waveform 
(|3.56|) is unique in that, unlike the formulas for evolving 
E, L z , it depends explicitly on this overall phase, and 
therefore on the initial position of the geodesic. 

At infinity, the effective energy-momentum tensor 
T^w is easily built from the GWs |gl|. The result can 
be expressed as [60| 



1 af3 



1 

lOTT 



dh+ dh A 



dh x dh y 



dx a dx$ dx a dx? 



(3.57) 



Brackets indicate an average over several wavelengths of 
the gravitational waves. Using this result and the expan- 
sion H3.56fl , it is straightforward to show that waves carry 
an averaged flux of energy and angular momentum given 

by 



dE 
~dt 



= Y — 

,u AlTUJ mkn 
Imkn mhii 

dL z \°° _ sr-^ 
~df I ~ a„,., 3 



lz H I 

2 | Imkn \ ' 



i mkn ^mkn 



J lmkn 



(3.58a) 
(3.58b) 



Similar expressions can be found for the flux of energy 
and angular momentum through the horizon; we refer 
the reader to Refs. Q, |5^, E2 for a detailed derivation. 
The resulting fluxes are 



dE\ H _ ^ 1 

dt J ^ 4tt/,i- 



imkn 



/^^p, a lmkn \ Zj lmkr. 



(3.59a) 



m kii 



dl^ 
dt 



47TW 3 ■ 

Imkn mKn 



The superscripts we use here are somewhat confusing: 
the fluxes to infinity depend upon the horizon coeffi- 
cients, Z^jfc n , and the fluxes down the horizon depend 
upon the infinity coefficients Z^ kn . This seemingly ob- 
tuse convention comes from the Green's function solution 
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we constructed, Eq. i|3.45|l . The coefficients Z n set the 
amplitude of the radial behavior towards infinity; the co- 
efficients Z°° set the amplitude towards the horizon. 
The coefficient ai m kn appearing in Eq. 1)3. 59[) is given 

by i 

_ 256(2Mr + ) 5 P(P 2 + 4e 2 )(P 2 + 16e 2 )^ fcra 

&lmkn " \f~f 10 7 

| (-' Imkn | 

(3.60) 

with e = v 'M 2 — a 2 /4Air+, where r+ is the location of 
the event horizon, and 

\Ci m kn\ 2 = [(A; m fe„ + 2) 2 + Aau m kn — 4a 2 w mfen ] 
x ( A Lfe« + 36maw mfe „ - 36a 2 w^ fc „) 
+ (2A /mfc „ + 3) (96a 2 cJ^ fcn - A%mau) mkn ) 
+ 144 W 2 nfe „(M 2 -a 2 ). (3.61) 

Recall that P = u) m kn ~ ma/(2Mr+), and that Xi m kn is 
related to angular equation's eigenvalue [cf. Eq. (|3.5b() ] 
via Eq. H3.14f) with lo = u) m kn ■ 

To go from snapshot waveforms to adiabatic inspiral 
waveforms, we must evolve not just energy and angular 
momentum, but also the Carter constant Q. This will 
eventually be done with the formalism recently derived 
by Mino [5J. However, if we assume that radiation does 
not influence an orbit's inclination [l or equivalently #i nc ], 
we can can infer an expression for (dQ/dt). This scheme 
was first suggested by Curt Cutler [63j]. From the def- 
inition (|2.7|) of l, we see that setting (dt/dt) = gives 

/dQ\ 2Q/dL z \ 

\* v = j:\-dt ■)' (3 - 62) 

where (dL z /dt) = {dL z /dt) H + {dL z /dt)°° . Hughes 
showed that this approximation is reasonable in the case 
of circular orbits |36j. It is currently unknown how ac- 
curate it may be in the case of generic orbits (though 
it is almost certainly pathological for orbits with nearly 
vanishing L z |6J|). In Sec. IVl below, we report values of 
(dQ/dt) given by Eq. (|3.62l) . We should emphasize how- 
ever that, since Eq. I|3.62|l is itself a speculation, those 
results will very likely be less accurate than the associ- 
ated fluxes of energy and angular momentum. 

IV. NUMERICAL ALGORITHM 

For a given orbit, the calculation of the snapshot wave- 
form is essentially identical to the calculation of each of 
the radiative fluxes. Given an engine which can com- 
pute the Z* mkn coefficients, which we will refer to simply 
as modes, each calculation is essentially just a long four 
dimensional summation. In the following three subsec- 
tions, we describe the details of our algorithms for (i) 
computing the frequency domain representations of the 
geodesies, (ii) computing the individual modes Z* mkn , 
and (iii) truncating the sums representing the waveforms 
and fluxes. 



A. Geodesies 

We begin by computing the quantities associated with 
the specified geodesic for a given set of black hole and or- 
bital parameters (a, e,p, 9i nc )- This calculation proceeds 
as follows: 

1. Compute the constants of the motion E, L z , Q 

2. Compute the orbital frequencies T r e,0, and T 

3. Construct routines to evaluate w r (ip), wg(x) and 
their derivatives 

4. Construct the expansions l|2.23a|) of At and A<fi 

Explicit expressions for the first three of these steps are 
shown in Appendix ^ The last step requires that we 
approximate the sums for the expansions (|2.23a() of Ax = 

At, Acj) as 

K N 

Ax 6 Ax e k sin(fcwg), Ax r » Ax^ sin(nuv), 

k=l n=l 

(4.1) 

where we determine the values of K and N by requir- 
ing that 15 consecutive Axl' coefficients make fractional 
contributions to the sum which are less than some speci- 
fied accuracy; we typically find K,N < 50. Each of these 
calculations associated with the geodesies is done to a 
fractional accuracy of 10 -12 , and the results are made 
available to the remaining computations. The calcula- 
tions described so far are typically completed in a few 
seconds (using a single ~ 1 GHz processor). It's worth 
bearing in mind that these calculations need only be done 
once when computing the waveform and fluxes for a given 
orbit. 



B. Modes 

After computing geodesic quantities, we construct an 
engine for calculating the modes Z* mkn . This requires so- 
lutions to the angular equation, and to the homogeneous 
radial equation. (Our homogenenous solutions have an 
implicit coupling to the source, since they depend on 
the frequencies oj m km which are determined by the or- 
bit.) The routines for solving the homogeneous radial 
and angular equations were inherited from the code used 
in Ref. [J, where they are described in detail. We solve 
the angular equation to nearly machine accuracy, and 
we solve the radial equation to a fractional accuracy of 
min(10 -7 , £fl ux /100), where £fl ux is the overall fractional 
accuracy demanded of the waveforms and fluxes. 

When computing (|3.50|) , we treat the radial integral as 
the outermost integral [the reverse of how we have writ- 
ten Eq. (|3.50|l ] . After experimenting with a variety of nu- 
merical methods for evaluating these integrals, we found 
a Clenshaw-Curtis quadrature algorithm to be the most 
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FIG. 2: Modal energy flux through the horizon (superscript 
H, marked with x's) and at radial infinity (superscript oo, 
marked with +'s). Results are for I = m = 4, k = 0. The 
orbit has eccentricity e = 0.7, inclination #i nc = 45°, and 
semilatus rectum p = 5.1, while the massive black hole has 
a = 0.9M. 



efficient. This method analytically integrates a numerical 
series representation (in Chebyshev polynomials) of the 
integrand, to produce a series representation of the inte- 
gral (for more details, see for example Ref. |65j). Each of 
these calculations depends on the indices (/, m, fc, n), and 
must be repeated many times to compute the waveform 
and fluxes associated with a given orbit. 

We compute the angular integral in Eq. I|3.50|) to a 
fractional accuracy of min(10 -6 , £fl ux /10). The accu- 
racy demanded of the radial integral is chosen dynam- 
ically. We begin by asking for a fractional accuracy of 
min(10 -5 , £flux). As each mode is computed, we store the 
magnitude of the largest modes encountered so far. Later 
modes are then computed to a fractional accuracy of at 
least 10%. After this 10% accuracy has been achieved the 
integrator continues to add more terms to the Chebyshev 
expansion. At each iteration it estimates the associated 
uncertainty in the total energy fluxes (the energy flux at 
the horizon for the Z^ kn integrals, and the energy flux 
at infinity for the Zf mkn integrals) by comparing to the 
magnitudes of the largest known modes. It then trun- 
cates the Chebyshev expansion as soon as the integral's 
associated relative flux error is smaller than £fl ux . The 
reason for this dynamic accuracy control is that many of 
the modes are very near zero. If we know that the total 
flux contains modes which are ~ 10 -3 , there is no need to 
compute a mode which is ~ 10 -20 beyond the first couple 
of digits. A single ~ 1 GHz processor usually takes con- 
siderably less than a second to compute both Z^ kn and 
Zfrnkn ( see Fig- 1 in Ref. |||), unless the modes are espe- 
cially badly behaved in the sense described in Sec. II V CI 
below. In these rare cases the modes can take anywhere 
from 10 seconds to a minute to compute. 



When computing a flux, the majority of the numerical 
work is devoted to computing modes which are insignif- 
icant or zero. The integrator has little difficulty when 
computing the dominant modes. However, modes which 
are very near zero (near round off error) can cause the 
integrator to fail. This typically happens when comput- 
ing modes with "large" values of the indices (l,m, k, n); 
the meaning of "large" is somewhat orbit and index de- 
pendent. For example, for fixed (l,m,k), we find that 
there is some limiting value of n beyond which the mode 
calculations cannot be trusted. 

We show an example of this behavior in Fig. El Beyond 
the dominant modes at n ~ 20, the mode magnitudes 
fall (roughly) exponentially until reaching a bottom. For 
the case shown in Fig. [21 the horizon modes bottom out 
at about n — 70, after which they gradually diverge, 
while the infinity modes bottom out at about n — 90, 
after which they appear to be dominated by numerical 
error. Both behaviors are unphysical since they indicate a 
divergent total energy flux. The different behaviors arise 
from the different scalings of these modes with frequency: 

\~JZ) ~ u mkn\ Z hnkn\ 2 > ( 4 - 2 ) 

\ Ui I Imkn 

\~jT/ ~ ^mkn\ Z WakS ' (4-3) 

\ 0,1 I Imkn 

This difference comes from the factor ai m k n ~ ^mkn m 
Eq. 

For reasonable requested accuracies, the modes which 
are difficult to compute are always insignificant when 
compared to the dominant modes (as is clearly the case 
in Fig. |2J). In the vast majority of cases, the truncation 
rules (see Sec. lIVCl belowl have been invoked long before 
these badly behaved modes are encountered. If this is not 
the case, these badly behaved modes can undermine the 
truncation scheme since, despite their small magnitudes, 
they do not regularly decay with increasing index. In or- 
der to account for this, we artificially zero any modes for 
which the ratio of the magnitude of the integral Zf mkn 
to the maximum over £ and x of the absolute value of 
the integrand in Eq. (|3.50() is below some threshold. We 
have found that this ratio is a good indicator of when the 
modes have bottomed out in the sense shown in Fig. [21 
Unfortunately, we have had to "hand-tune" this param- 
eter somewhat. This was particularly so as we explored 
orbits with increasing eccentricity. In the future, we ex- 
pect to incorporate our code into one which computes 
inspiral waveforms. Hand tuning of this parameter will 
be unacceptable for such a code, and a more robust tech- 
nique will be needed. 

C. Truncation 

Given the routines described above, the remaining task 
is to evaluate the four dimensional sums which make up 
the waveforms and fluxes (Sec. IIIID|I . As mentioned 
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above, we do this by monitoring the sums which rep- 
resent the fluxes. We assume that once the fluxes have 
converged to some requested fractional accuracy £fl U x, the 
sums representing the waveforms are comparably accu- 
rate since they are constructed from the same mode co- 
efficients as were used for the fluxes. We now discuss the 
algorithm for truncating the sums which represent the 
fluxes. 

Let F be the flux of energy or angular momentum at 
the horizon or at infinity 



H,oo 



dL_2 
dt 



H,oo 



The algorithm for computing these fluxes is: 

oo 



1=2 
I 



Fi= Flm > 



Fl m — 2_j ^ 



link : 



k= — oo 



Flmk — F[ 



mkO 



2j> 



rnkn • 



(4.4) 

(4.5a) 
(4.5b) 
(4.5c) 
(4.5d) 

(4.6) 



We have used the symmetry Ij3.52l) . which gives 

Flmkn — -F 1 ;(-m)(-fc)(-n); 

to simplify the sum over the radial harmonic index n. 
Each of the sums which has an infinite boundary is then 
truncated as described below. 

The terms which make up outermost sum, F[ in the 
sum over I, always decrease monotonically with increas- 
ing I. This means that the Z-sum can be written 



F = J2 F i+ 6F > 



(4.7) 



i=i 



where to a good approximation SF w Fl, and I is a 
number which is determined by the requested accuracy 
via 



\F L \ < £fi ux max|F ; | 



(4.8) 



The sum over the m-index is computed for all values 
-I < m < I — we do not truncate it. 
The k and n sums are approximated as 



K 4 



JV 4 



F[m ~ ^ ^ Fi 7n kj -^Zmfc ~ Fi 7n kQ -\- 2 ^ ^ -^lm/cm 



k=K_ 



n=l 



(4.9) 



where the constants K± and N + are determined by re- 
quiring that the following relation 

\FimK ± \ < £flux max l-FKn'fe'l, (4.10a) 

I'm'k' 

\FimkN+\ < £flux max |F;/ 

m' k'n' 

|, (4.10b) 

I'm'k'n' 



(where the maximization is done over all previously en- 
countered mode amplitudes) is satisfied for B consecutive 
terms. We also require that the first of these B terms is 
larger in magnitude than the last. This extra condition 
gives some assurance that we are not truncating during a 
slow rising up of the modes. For the sums over the polar 
harmonic index fc, we use B = 2. For the sums over the 
radial harmonic index n, we use B = 5 — the dominant 
radial modes show a less orderly distribution than the 
polar modes. 

The validity of any given truncation scheme is always 
subject to doubt — it could be the case that modes which 
haven't been computed are anomalously large in magni- 
tude. In this sense truncating sums can be likened to 
financial investments: past performance is not necessar- 
ily indicative of future results. 

We have found these words of caution are especially 
applicable for the sums over the radial harmonic index 
n. We find that the Fi terms, as defined by Eqs. (|4.5|l . 
in the outermost sum over I always decrease monotoni- 
cally with increasing I. Also, the Fi m k terms are always 
peaked near k — 0, so that setting the buffer Btoa mod- 
estly small number is sufficient for catching the dominant 
modes. Unfortunately, the Fi m k n appear not to obey any 
similar general rules. If either I or the eccentricity e are 
increased, the n values of the dominant Fi m k n modes drift 
to large values of n. As a result, the truncation rules are 
expected to be limited in a rough sense by large eccen- 
tricity, or by requests for high accuracy (which would 
require modes with large I). 

It may be possible to better anticipate n for the 
dominant modes. Peters and Mathews [|3(j used the 
quadrupole approximation to describe the radiation of 
two point particles moving along closed Keplerian orbits 
of arbitrary eccentricity e. Under this approximation, 
they decomposed the metric perturbation into tensor 
spherical harmonics [67| with 1 = 2. They then showed 
that the power (diJ/dt)™ radiated at a frequency nio 
(where n is an integer and to is the single unique orbital 
frequency for a Keplerian orbit) is given by 

PM 4 



/dEY m n 4 f fT , , , . 

\~dt / K 32 \ l Jn - 2 ( ne > ~ 2eJ n-i(ne) 



+ -J n (ne) + 2eJ n+ i(ne) - J n+2 (ne)] 
n 

+ (1 - e 2 ) [J„_ 2 (ne) - 2J n (ne) + J n+2 {ne)] 2 

4 2 
'■in 



(4.11) 



where J n (x) are Bessel functions of the first kind. We 
have found that this formula (which was derived using 
only the quadrupole equation and Keplerian orbits) accu- 
rately predicts the n values of the dominant Fi m k n modes 
for 1 = 2. More quantitatively, define h by 



dE\ PM I dE \ PM 

— — ) = max ( — — ) 

«/fi " \dt/ n 



(4.12) 
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FIG. 3: Modal energy flux through the horizon (superscript 
H, marked with x 's) and at infinity (superscript oo, marked 
with +'s). The indices I = m — 2, k = 0. The orbit has 
eccentricity e = 0.85, inclination lnc = 20°, and semilatus 
rectum p = 6, while the massive black hole has a = 0.9A/. 
The horizon modes are peaked at n — 29, while the infinity 
modes are peaked at n = 26. The dashed line shows the 
estimated peak location (at n — 28) found using the Peters- 
Mathews power formula (|4. lipi . 



so that 



n « exp 



(4.13) 



this approximate formula was found by fitting to numer- 
ically determined maxima of Eq. (|4.11l) . Then, even for 
highly relativistic, inclined, and eccentric orbits, we have 



max \F- 



2mkn \ 



\F. 



2mkh | 



(4.14) 



An example is shown in Fig. EI Reproducing Fig. [3] un- 
der the same conditions but with an orbital inclination 
of #inc = 70° shifts the peaks occur to n ~ 35. It is some- 
what remarkable that the Peters-Mathews result can be 
used to predict the peak locations in this way, since the 
generic Kerr orbits considered here do not closely resem- 
ble the analogous Keplerian orbits either in their orbital 
frequencies or in their shapes. We note that while the 
Peters-Mathews power formula (|4.11|l gives reasonable 
estimates for the location in n of the dominant modes, 
it does not accurately predict the actual flux associated 
with the dominant modes (for example, the top panel of 
Fig. 14 in Ref. shows a very narrow peak when com- 
pared to the relatively broad peak for the e = 0.7 orbit 
in Fig. 3 from Ref. [66^ . 

It is likely that the I > 2 analog of the Peters-Mathews 
power formula (|4.11() would be a useful estimator of the 
values of the radial harmonic index n for the dominant 
modes. Such an estimator would be valuable tool for 
designing a more accurate truncation scheme for the sums 



over n. The value of n for the dominant modes Fi m k n 
increases roughly linearly with I, so it may be possible 
to "track" the peak location starting from the Peters- 
Mathews power formula 14. lip at I = 2. Future progress 
along these lines may reduce the limitations of our code 
for computing accurate fluxes for highly eccentric orbits. 



D. Validation 

We have compared the output of our code with two 
existing codes. The first, the "circular code" by Hughes 
0, treats orbits which have a non-zero inclination and 
a constant radius. The second, the "equatorial code" 
by Glampedakis and Kennefick 0, treats orbits which 
are eccentric and confined to the equatorial plane. Both 
of these codes have already passed a variety of tests of 
their own (described below, but see Refs. 0, Q for more 
details). Also, the circular and equatorial codes have 
been shown to agree with each other for orbits which are 
both circular and equatorial 0|. 

The circular code has been shown to agree with analyt- 
ical post-Newtonian approximations, found in Ref. |59| . 
in the weak field limit. It has also been shown to have 
the correct rotational properties in the Schwarzschild 
limit (see Sec. V B of Ref. Ij for details). Our code 
is a direct descendant of the circular code, and it has 
some elements in common with that code. Namely, both 
use the same routines for solving the homogeneous ra- 
dial equation, and also for solving the angular equa- 
tion. However, they differ in their methods for solving 
the geodesic equations, in their methods for computing 
the quantities Zf mkn , and in their methods for truncat- 
ing the sums representing the fluxes 

\H,oo 



When comput- 
ing the individual modes {dX/dt)^°^ (for X = E,L Z ) 
with our code and the circular code, we find differences 
A {dX/dt)f m ™ n < (1(T 5 ) (dX/dt) U -°°. Here we are com- 
paring the difference to the total flux rather than the 
mode itself because often there can be large fractional 
disagreements for insignificant modes. When computing 
the total fluxes (dX/dt} ti ' co with these two codes, we find 
differences A (dX/dt)^ 00 < (1CT 5 ) {dX/dt) n '°°. 

The equatorial code has been tested against its own 
direct ancestors [U|68|, against the circular code, and 
also against Shibata's earlier equatorial code [3^. Note 
that in the case of the test against Shibata's code, 
Glampedakis and Kennefick found agreement only at a 
level of about 1%. Shibata has since reported that after 
improving his code, he now agrees with Glampedakis and 
Kennefick to much greater accuracy |69| . When comput- 
ing the individual modes (dX/dt) lr ^ n (for X = E,L Z ) 
our code always agrees with the equatorial code at the 
level of the equatorial code's claimed fractional accu- 
racy of 10 -3 . We have tested this for the top panels of 
Figs. 14-16 of Ref. 0; by fractional accuracy, we mean 
the ratio of the differences in modes to the sum of all the 
modes shown in each of the figures. When computing the 
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total fluxes (dX/dt) we find that our code always agrees 
with the equatorial code at the level of their claimed frac- 
tional accuracy of 10 -3 . We have tested this for the top 
panel, for a = 0.5, of Table VII of Ref. 0. Note that 
the fractional accuracy of 10 -3 is the most conservative 
of the accuracy claims made in Ref. 0. Usually we find 
agreement with the equatorial code well beyond this fig- 
ure; for many of the total fluxes, we agree with all six of 
their published digits. 

While we find agreement with the equatorial code when 
computing the total fluxes (dX/dt), the two codes show 
significant differences when computing the horizon fluxes 
(dX/dt) 11 (no such disagreement was found for the fluxes 
at infinity (dX/dt)°°). The disagreement is especially 
strong for the energy flux at the horizon, but it is also 
present to a lesser extent with the angular momentum 
flux at the horizon. For the energy fluxes at the horizon, 
we typically only agree with the equatorial code to about 
1-10%, however in one case (the orbit with e = 0.4 and 
p = 5 in Table VII of Ref. 0] , which has an exceptionally 
weak energy flux at the horizon) our results differ frac- 
tionally by about a factor of three, and they even have 
different signs. Currently the source of this disagreement 
is unknown. We emphasize however that the total flux of 
either energy of angular momentum is always dominated 
by the flux at infinity, and even in the most extreme case 
of our disagreements over the horizon fluxes, we still find 
acceptable agreement for the total fluxes. 

Hughes has shown that, in the Schwarzschild limit, un- 
der a transformation from an equatorial to an inclined 
orbit, the individual modes must transform as |l| 12 



{dE/dt)l 



l{m-k)kn K v mc 



(dE/dt); m0n (e 



o) 



V 



(m-fc)™-^ inc 



(4.15) 



where T3 l m , m {0 mc ) is the Wigner D-function, as defined 
by Eq. (4.256) in Ref. |H: 



min(Z— m,l-\-m f ) 



X COS ■ 



q— max(0, m' —m) 
2l-\-7n' —m—2q 



m — m -\-2q 



+ m')\{l - mfjUJ + m Y-( l - TO ) ! 
q\{l — m — q)\(l + ml — q)\(m — ml + q)\ 



(4.16) 



For a variety of generic Schwarzschild orbits (which are 
both eccentric and inclined) and for a variety of modes, 
we have checked that our code satisfies this require- 
ment to a fractional accuracy of 10~ 5 or better. We 
have also confirmed that, again for a variety of generic 
Schwarzschild orbits, the total fluxes of energy at both 
the horizon and infinity, are independent of inclination 



12 Here we correct two typos in Eq. (5.4) of Ref. Qj. 
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FIG. 4: A comparison of waveforms associated with a generic 
orbit and with the related circular and equatorial orbits. Each 
orbit has a semilatus rectum of p — 4. The magnitude of the 
massive black hole's spin angular momentum is aM = 0.9Af 2 . 
Each waveform is plotted for 1.5 x M / (W 6 Mq) hours as seen 
from a viewing angle of 9 = 30° away from the symmetry axis 
of the massive black hole. 



(we have checked this up to a fractional accuracy of 

io- 4 ). 

With the help of Norichika Sago, we have also 
compared the output of our code with analytical 
post-Newtonian expressions (essentially an extension of 
Ref. [5^) for the fluxes of energy and angular momentum 
at infinity for orbits which are both slightly eccentric and 
slightly inclined 70] . We compared the contributions for 
I = 2, 3, and 4. With a requested accuracy of 10~ 4 
in our numerical code, the two methods of calculation 
always agreed to ~ 10~ 3 or better (with slightly bet- 
ter agreement for angular momentum fluxes than for en- 
ergy fluxes). Preliminary investigations suggest that this 
disagreement is due to the small-eccentricity-expansion 
used when deriving the post-Newtonian expressions, and 
is therefore not unexpected. 



V. RESULTS 

In this section we discuss the results found from using 
the code described in the previous section. After de- 
scribing the general characteristics of the radiation from 
a typical generic orbit, we will discuss how these charac- 
teristics vary over a catalog of orbits. 



A. Radial and polar voices 

Figure ^ compares a snapshot waveform produced by 
a generic orbit with the waveforms produced by the cir- 
cular and equatorial limits of that orbit. It is clear that 
the generic waveform is not well approximated by either 
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of the limiting cases at the level of accuracy needed to 
produce detection waveforms (fractional phase accuracy 
of ~ 10~ 4 ). We can however understand the generic 
waveform as a composition of the effects produced by 
the different components of the orbital motion: the polar 
motion 9(X) and the radial motion r(A). To do so, we 
first define the general waveform quantity 



H = h + - ih x = 2J H kn t 



-iu) kn {t~r) 



kn 



where 

Hkn = _ ~E 



9 7 H 

mhn 



(5.1) 

(5.2) 
(5.3) 



Each term Hk n exp[— iu>k n {t — r)] in the sum l|5.1[) can 
be likened to a "voice" in the GW produced by the or- 
bit. By using the general expression for the Z* mkn co- 
efficients <|3.5(J|1 and the symmetry rules (I2.24[) for the 
geodesic expansions, we see that in the limiting circular 
and equatorial cases, only some voices contribute 

Hkn oc S n o, (circular orbits) (5.4a) 
Hkn oc 5ko, (equatorial orbits) (5.4b) 
Hkn SnoSko- (circular-equatorial orbits) (5.4c) 

The term -ffoo, the sum of all the terms which oscil- 
late at pure multiples of the azimuthal frequency, is the 
only nonvanishing term for orbits which are both circu- 
lar and equatorial. The polar and radial voices are made 
of terms which oscillate at integer linear combinations of 
the azimuthal frequency, and either the radial or polar 
frequency. The polar voices, with n = and k ^ 0, can 
be thought of as being produced by the orbital inclina- 
tion. Similarly, the radial voices (k — and n / 0) are 
produced by the orbit's eccentricity. This suggests the 
following separation of the voices: 

H ^azimuthal ~t~ ffpolar -\- H radial "h -ffmixcdi 

(5.5) 

where 



-ffazimuthal 
-ffpolar 

-ffradiai 

-ffmixed 
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fc#0 



n=£0 



HonC 



— iujQ n t 



EE^ 



(5.6a) 
(5.6b) 

(5.6c) 

(5.6d) 



Figure displays the waveforms for the different voices 
of the generic waveform from Fig. 0] Since that orbit 
has eccentricity e = 0.5 and inclination 9[ nc = 45°, one 
might expect that the radial and polar voices would be 
comparably loud. Instead (as one might guess by looking 
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FIG. 5: The waveform from Fig. 2] and its separated pieces 
Eq. I|5.6a|l . This orbit has an eccentricity e = 0.5, inclination 
#inc = 45°, and semilatus rectum p = 4. The spin of the 
black hole is a — 0.9M. We plot the waveforms as seen from a 
viewing angle of 9 = 30° . A plot of h x shows similar behavior. 





radial 


polar 


azimuthal 


mixed 


(dE/dt) 


62% 


1.2% 


0.52% 


37% 


(dL z /dt) 


73% 


0.77% 


1.0% 


25% 



TABLE II: The results of applying the splitting procedure 
described by Eq. 15.51 to the fluxes of energy E and angular 
momentum L z produced by the generic orbit from Fig. [1] 
Each number is the ratio of the contribution from a specific 
voice, to the total quantity. 



at the waveform from the equatorial limit of this orbit; cf. 
Fig. the radial voice alone -ffradiai has a phase history 
which is similar to that of the complete waveform. The 
polar voice H po i m - is comparatively unimportant. 

The fluxes in energy and angular momentum can be 
similarly separated into radial, polar, azimuthal, and 
mixed pieces. Again for this case, the radial voice dom- 
inates. Table ITT1 shows how the fluxes are distributed in 
the case of the orbit used to produce Figs. 0] and The 
radial voice carries away more than half of both the en- 
ergy and angular momentum fluxes, while the polar voice 
carries only about 1%. This suggests that GWs from this 
orbit might, for some applications, be well approximated 
by its radial voice alone: H « -ffradiai- Such an approxi- 
mation would by highly desirable since it would consid- 
erably reduce the computational cost of generating this 
waveform. Unfortunately -ffradiai is still quite sensitive to 
the orbital inclination, as can be seen by comparing to the 
waveforms for the circular and equatorial orbits in Fig. 0] 
This is because the inclination influences the values of all 
three orbital frequencies. So even if this approximation 
were valid for some range of orbits, it would not likely 
reduce the number of detection waveforms needed when 
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searching for events produced by those orbits, though 
it may reduce the computational cost to produce each 
waveform. 



B. Catalog of orbits 

We now discuss the waveforms and fluxes for a cata- 
log of generic orbits. The parameters of these orbits are 
shown in Table II I II Tables IIVI and [V] show the relative 
contributions of the different voices for these waveforms. 

It is clear that as inclination is increased, radiation 
is channeled into the polar voice, and that likewise the 
radial voice is amplified by an increase in eccentricity. 
Although both voices can be amplified by such adjust- 
ments, the radial voice seems particularly booming while 
the polar voice is more subtle. While the polar voice can 
be made dominant, it requires especially low eccentric- 
ity, below about 0.3. Also, for orbits with eccentricity 
as low as 0.1 and inclinations mc < 45° or 9 lnc > 135°, 
one should expect to capture at least half of the radia- 
tive power and torque in the azimuthal voice alone. This 
is especially significant since the azimuthal voice is ex- 
ceedingly inexpensive to compute in comparison to the 
others. It is composed only of oscillations at pure multi- 
ples of the azimuthal frequency with k = n = 0, so that 
computing it alone would mean replacing all of the four 
dimensional sums in Sec. lIII Dl with only two dimensional 
sums. 

The asymptotic energy and angular momentum fluxes 
associated with the catalog ( Table ITTTf) are displayed in 
Tables EH and ETD The energy and angular momen- 
tum fluxes were requested to have fractional accuracies 
of £flux = 10~ 4 . The estimates of the actual fractional 
accuracies, shown in square brackets in Tables IVII and 
IVIII were determined as follows: First we compute each 
flux with a requested fractional accuracy of £fi ux = 10 -4 . 
The fractional accuracy claimed in Tables IVII and IVIII is 
then either 10~ 4 , or if larger, the fractional residual that 
was found when comparing with the same flux computed 
with a requested fractional accuracy of £fl ux = 10 -3 . The 
summation buffer B for the polar harmonic index k is 3, 
while B — 5 for the radial harmonic index n. The largest 
I value needed for the e = 0.1, 0.3 orbits was I — 10, while 
for the e = 0.5, 0.7 orbits the largest I value was I = 11. 
The most complex waveform, (e,6*i nc ) = (0.7,80°), is 
made up of about 160,000 modes (about 80,000 mode 
calculations due to symmetry). The simplest waveform, 
(e,0inc) = (0.1,20°), is made from about 17,000 modes 
(about 8,500 mode calculations). 

Some examples of waveforms are finally shown in the 
time domain in Figs. El There we plot h + waveforms for 
3 x M / (10 6 Mq) hours as seen from four different viewing 
angles 6> = 0°, 30°, 60°, 90°. These plots show that the 
radial voice has a burst-like character, while the polar 
voice is more of a modulated hum. These features are 
consistent with the characteristics found in earlier work 



which focused on circular orbits jl| and on equatorial 
orbits 0. 



VI. SUMMARY AND FUTURE WORK 

This work describes the first calculation of gravita- 
tional waves, and asymptotic radiative fluxes of energy 
and axial angular momentum, produced by a spinless test 
mass on a generic bound geodesic of a spinning black 
hole. It represents the natural extension of the earlier 
works shown in Table |U The future direction of this 
work is rather straightforward: fill in the last checkmark 
in Tabled by implementing an evolution scheme for the 
constants of the geodesic orbits, and in particular for the 
Carter constant. Our anticipatedpath toward this goal is 
discussed in more detail in Ref. [g. Once completed, we 
will be able to compute generic inspiral waveforms which 
will facilitate the initial detection of EMRIs. Completing 
this program will require implementing Mino's scheme for 
evolving the Carter constant into equations which can be 
entered directly into a code H, IE S 13 ■ ^ wu ^ a l so re ~ 
quire a variety of more computational tasks, such as de- 
veloping a dynamical scheme for artificially zeroing badly 
behaved modes as described in Sees. IIV O and llV Bl We 
are optimistic about completing this program since most 
of the remaining work is relatively straightforward in the 
sense that there are no known potentially overwhelming 
hurdles. 
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TABLE III: Each row above describes an orbit in our catalog. The numbers to the right hand side of the vertical bar were 
computed (to an fractional accuracy of 10 -12 ) from the numbers on the numbers on the left hand of the bar (using the method 
described in Appendix lA^ . Each orbit has a — 0.9M. 



through the internal Research and Technology Develop- APPENDIX A: GEODESIC DETAILS 

ment program. 

The functions Vt, V r , V$, and V$ appearing in the 
geodesic equations (|2.4|) are given by 

V t (r,6) 

Vr(r) 
Vg{9) 



= E (x ~ 0,2 sin2 e ) + aLz i 1 ~ x) ' 

(Ala) 

= {Em 2 - aL z f - A [ M V + [L z - aE) 2 + Q] , 

(Alb) 

= Q-L 2 z cot 2 6-a 2 (fi 2 -E 2 )cos 2 6, (Ale) 
= L z esc 2 9 + aE f^- - l) - (Aid) 
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TABLE IV: Here we show the results of applying the splitting procedure described by Eq. 1)5. 5fl to the fluxes of energy E and 
angular momentum L z produced by our catalog of generic orbits (Table IrTTt . The symbol V represents a ratio of the average 
power relative to the total average power {dE/dt), and the symbol T represents a ratio of the average torque to the total 
average torque (dL z /dt). Each of these orbits has a — 0.9M and p — 6. 
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TABLE V: Identical to Table ITvl except that these orbits are retrograde and have semilatus rectum p = 12. 



where we have defined 

m 2 =r 2 + a 2 . (A2) 

We now show some explicit relations which are used to 
evaluate functions of the geodesies. As indicated below, 
these relations were first derived either by Schmidt 
or in Ref. j^. We reproduce them here for the sake of 
completeness. 

When evaluated in order, the following equations de- 



termine the energy E = fiE, angular momentum L z = 
fiML z , and Carter constant Q = (fiM) 2 Q for an orbit of 
a given minimum radius f\ — r m ; n /M, maximum radius 
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TABLE VI: The fluxes of energy E and axial angular momentum L z at infinity (superscript of oo) and through the black 
hole's event horizon (superscript of H) for some generic orbits. Each of these orbits has a = 0.9M and p = 6. Each number in 
square brackets is an order of magnitude estimate for the fractional accuracy of the preceding number. The rates of change for 
the Carter constant (dQ/dt) were computed from the formula 13.621 which follows from the assumption that radiation does not 
change the orbital inclination. Since this is an uncontrolled approximation, the accuracy of the (dQ /dt) figures is unknown; we 
include these data for comparison purposes. 
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TABLE VII: Identical to Table IvTl except that these orbits are retrograde and have semilatus rectum p — 12. 



^2 = r ma x/M, and orbital inclination #i llc = tt/2 — 9 m - m : These were derived by Schmidt in Appendix B of 

Ref. Here the constant D — sgnL z , so that D = 1 

~ 2 np + 2ea— 2Dy/ a(cre 2 + pen — rjn 2 ) for a prograde orbit, and D = — 1 for a retrograde or- 

J " -2 i i \ ) bit. The other constants in these equations are defined 
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FIG. 6: Snapshot waveforms for orbits with inclination 8i nc = 80°, semilatus rectum p — 6, and eccentricities e = 0.1, 0.3, 0.5, 
0.7. The magnitude of the black hole's spin angular momentum is aM = 0.9M 2 . 



in terms of the constants 13 

a = a/M, z_ = cos 2 min , (3 = a 2 (l - E 2 ), (A6) 
the functions 



and the determinants 
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(A13) 
(A14) 
(A15) 
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A subscript 1, 2 means the function is to be evaluated at 
ri,2- 



Note that Schmidt defines z_ as cos 6* m ; n l45l . 



Schmidt also derives the following expressions for the 
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coordinate-time frequencies: 



where the functions F, G, H, J are given by: 



MCLr = 



■npK(k) 
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Here k = \J z-/z+, 



z + 



L 2 Z + Q + f3 + \f (L 2 Z + Q + /3) 2 - 4f3Q 
~2[3 ' 



(A20) 



K(k) is the complete elliptic integral of the first kind, 
E{k) is the complete elliptic integral of the second kind, 
and n[<£, n, k) is the Legendre elliptic integral of the third 
kind 65 1 : 



K{k) = 
E(k) = 
H(4>, n, k) = 



tt/2 



(10 
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o (1 - n sin 2 9) \/ 1 - k 2 sin 2 i 



(A21) 
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(A23) 



The remaining quantities in Eqs. l)A17jl - lA19|) are defined 
by the following integrals: 



W = 
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Y = 
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(A31) 

The Mino time frequencies are then found using T 0,0, r = 
1^,^, and T e = tt v7^+/[2if(fc)] (from Ref. Qj. 

The function uy (•;/> ) a nd its derivative were derived in 
the Appendix of Ref. J44j; however, there we did not write 
things explicitly in terms of Schmidt's J- function. The 
results are: 
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The function uuni x) a nd its derivative was derived in the 
Appendix of Ref. [iij . The results are 



we{x) 
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2K{k) 1 - fc 2 cos 2 x' 
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(A36) 



and T(<b, k) is the incomplete elliptic integral of the first 
kind 1691: 



H4>,k)= I -= 

Jo VI 
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\Jl - k 2 sin 2 



(A37) 



20 

APPENDIX B: PERTURBATION DETAILS The A abc functions are given by 14 



The differential operators in the master equation (|3.2() 
are given by 
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(Bib) 



In Boyer-Lindquist coordinates, the Kinnersly tetrad 
vectors are given by Q, 
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and the Kinnersly tetrad one-forms are 
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A ~ 1 ~ aAsin 2 
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14 Note that the code used in Refs. Hl l ifl had the wrong prefactor 
for A nn o (it had the factor of 1/p replaced with a factor of 1/p). 
Also Ref. Q has an incorrect expression for A n mi (though this 



where we have used the shorthand notation S = 
Si m (9,uj). Note that the factor of C^m in Eq. (|B4c|) 
was missing in Ref. pj (it was not however missing from 
the numerical code associated with Ref. Q). 

term is correct in the code for Refs. HI l.'itl i 
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